;;----------------------------------------------------------------------
;; scatter plot 
;;----------------------------------------------------------------------
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "$KAILIBROOT/callvar.ncl" 
  load "$KAILIBROOT/callstr.ncl" 

  print("") 
  print("") 
  print("") 
  print("") 
  print("") 
  print("") 
;;----------------------------------------------------------------------
;; user setting 
;;----------------------------------------------------------------------

  colormap = "perc2_9lev" 
  colormap = "testcmap" 
  icla = 48  
  iclb = 60 
  iclc = 148 

  colormap = "psgcap" 
  icla = 148  
  iclb = 106  
  iclc = 109 
  iclc = 178 


  LeftString   = ") Radon (mBq m~S~-3~N~ STP)"  
  RightString  = ""  
  XAxisString  = "Observation ~C~  "
  YAxisString  = "Model"
  plotform     = "eps"  
  plotname     = "fig_scatter_color"  
  L_log = 0
  axisfont = 0.025 
  textfont = 0.018
  np = 9
  plot = new(np,graphic )

  markersize = 0.030 
  markertype = 1 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------
  ncfnama = "radon_mine_E0_1999-2003_ymm.nc" 
  ncfnamb = "radon_mine_E1_1999-2003_ymm.nc" 
  ncfnamc = "radon_mine_E2_1999-2003_ymm.nc" 
  ncfnamd = "radon_mine_E3_1999-2003_ymm.nc" 
  ncfname = "radon_obs_mm.nc" 
  ascnama = "radon_siteplot.dat" 
  ascnamb = "radon_sitename.dat" 
  ascnamc = "radon_sitebad.dat" 
  ascnamd = "radon_sitelat.dat" 
  ascname = "radon_sitelon.dat" 
  ascnamf = "radon_siteeur.dat" 
  ascnamg = "radon_sitechi.dat" 
  ascnamh = "radon_siteusa.dat" 
  ascnami = "radon_sitesam.dat" 

  fla = addfile(ncfnama,"r") 
  flb = addfile(ncfnamb,"r") 
  flc = addfile(ncfnamc,"r") 
  fld = addfile(ncfnamd,"r") 
  fle = addfile(ncfname,"r") 

  rma = fla->radon
  rmb = flb->radon
  rmc = flc->radon
  rmd = fld->radon
  rme = fle->radon

  kplot = asciiread(ascnama,-1,"integer") 
  sname = asciiread(ascnamb,-1,"string") 
  kbad  = asciiread(ascnamc,-1,"integer") 
  lat   = asciiread(ascnamd,-1,"float") 
  lon   = asciiread(ascname,-1,"float") 
  ieur  = asciiread(ascnamf,-1,"integer") 
  ichi  = asciiread(ascnamg,-1,"integer") 
  iusa  = asciiread(ascnamh,-1,"integer") 
  isam  = asciiread(ascnami,-1,"integer") 

  ns = dimsizes(kplot) 

  print("ns: "+ns) 

  do is = 0,ns-1
    if(ieur(is).eq.0) then 
       print("ieur: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(ichi(is).eq.0) then 
       print("ichi: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(isam(is).eq.0) then 
       print("isam: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(iusa(is).eq.0) then 
       print("iusa: "+sname(is)) 
    end if 
  end do 


  knot  = kbad 
  kval  = kbad 


  rma!1 = "site" 
  rmb!1 = "site" 
  rmc!1 = "site" 
  rmd!1 = "site" 

;;----------------------------------------------------------------------
;; mask out bad sites 
;;----------------------------------------------------------------------
  do is = 0,ns-1 
     if(kbad(is).eq.0) then 
        rma(:,is) = -999  
        rmb(:,is) = -999  
        rmc(:,is) = -999  
        rmd(:,is) = -999  
        rme(:,is) = -999  
     end if  
  end do 

  do is = 0,ns-1
     if(all(ismissing(rme(:,is)))) then
        knot(is) = 0 
     end if
  end do
 
;;----------------------------------------------------------------------
;; annual average 
;;----------------------------------------------------------------------
  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  

  nt = 12 

;;----------------------------------------------------------------------
;; mask out sites only have annual mean 
;;----------------------------------------------------------------------

  do is = 0,ns-1

     xrmaa = rma(0,is) + rma(1,is) + rma(11,is)  
     xrmab = rma(2,is) + rma(3,is) + rma(4,is)  
     xrmac = rma(5,is) + rma(6,is) + rma(7,is)  
     xrmad = rma(8,is) + rma(9,is) + rma(10,is)  

     xrmba = rmb(0,is) + rmb(1,is) + rmb(11,is)  
     xrmbb = rmb(2,is) + rmb(3,is) + rmb(4,is)  
     xrmbc = rmb(5,is) + rmb(6,is) + rmb(7,is)  
     xrmbd = rmb(8,is) + rmb(9,is) + rmb(10,is)  

     xrmca = rmc(0,is) + rmc(1,is) + rmc(11,is)  
     xrmcb = rmc(2,is) + rmc(3,is) + rmc(4,is)  
     xrmcc = rmc(5,is) + rmc(6,is) + rmc(7,is)  
     xrmcd = rmc(8,is) + rmc(9,is) + rmc(10,is)  

     xrmda = rmd(0,is) + rmd(1,is) + rmd(11,is)  
     xrmdb = rmd(2,is) + rmd(3,is) + rmd(4,is)  
     xrmdc = rmd(5,is) + rmd(6,is) + rmd(7,is)  
     xrmdd = rmd(8,is) + rmd(9,is) + rmd(10,is)  

  ;if(is.eq.0) then    ;; beijing
  ;   rme(:,is) =  -999 
  ;   rme(0,is) =  9078. 
  ;   rme(3,is) =  6248. 
  ;   rme(6,is) =  6847. 
  ;   rme(9,is) =  9945. 
  ;end if  

  if(is.eq.93 .or. \ 
     is.eq.94 .or. \ 
     is.eq.95 .or. \ 
     is.eq.96 .or. \ 
     is.eq.97 .or. \ 
     is.eq.98 .or. \ 
     is.eq.100 ) then 
     rma(:,is) =  -999
     rmb(:,is) =  -999
     rmc(:,is) =  -999
     rmd(:,is) =  -999
     rma(0,is) = xrmaa / 3.  
     rmb(0,is) = xrmba / 3. 
     rmc(0,is) = xrmca / 3. 
     rmd(0,is) = xrmda / 3. 
     rma(3,is) = xrmab / 3. 
     rmb(3,is) = xrmbb / 3. 
     rmc(3,is) = xrmcb / 3. 
     rmd(3,is) = xrmdb / 3. 
     rma(6,is) = xrmac / 3. 
     rmb(6,is) = xrmbc / 3. 
     rmc(6,is) = xrmcc / 3. 
     rmd(6,is) = xrmdc / 3. 
     rma(9,is) = xrmad / 3. 
     rmb(9,is) = xrmbd / 3. 
     rmc(9,is) = xrmcd / 3. 
     rmd(9,is) = xrmdd / 3. 
  end if 

  if(is.eq.93) then   ;; huhehaote
     rme(:,is) =  -999 
     rme(0,is) = 10600.
     rme(3,is) =  6900.
     rme(6,is) =-10000.
     rme(9,is) =  7340. 
  end if   
  if(is.eq.94) then    ;; changchun
     rme(:,is) =  -999 
    ;rme(0,is) = 10500.
     rme(0,is) = 13000.
     rme(3,is) =  8800.
     rme(6,is) = 11380.
     rme(9,is) =  9390. 
  end if  
  if(is.eq.95) then   ;; nanjing  
     rme(:,is) =  -999
     rme(0,is) =  6410.
     rme(3,is) =  6420.
     rme(6,is) =  6430.
     rme(9,is) =  6440.
  end if
  if(is.eq.96) then   ;; xi_an
     rme(:,is) =  -999 
     rme(0,is) = 17400. 
     rme(3,is) =  9600. 
     rme(6,is) =  9200. 
     rme(9,is) = 10000
  end if 
  if(is.eq.97) then   ;; shanghai 
     rme(:,is) =  -999
     rme(0,is) =  4810.
     rme(3,is) =  4820.
     rme(6,is) =  4830.
     rme(9,is) =  4840.
  end if
  if(is.eq.98) then    ;; wuhan
     rme(:,is) =  -999 
     rme(0,is) = 23920.
     rme(3,is) = 13770.
     rme(6,is) =  8500.
     rme(9,is) = 12670. 
  end if  
  if(is.eq.100) then   ;; guiyang
     rme(:,is) =  -999
     rme(0,is) = 12210.
     rme(3,is) = 12220.
     rme(6,is) = 12230.
     rme(9,is) = 12240.
  end if

  if(is.eq.30) then   ;; rio
     rme(1,is) = 1450.
     rme(8,is) = 1560.
     rme(9,is) = 1460.
     rme(11,is)= 1390.
     print(sname(is)+" "+rme(:,is)) 
  end if

  if(is.eq.29) then   ;; cha
     rme(1,is) = 1450.
     rme(3,is) = 1490.
     rme(4,is) = 1470.
     rme(5,is) = 1600.
     rme(7,is) = 1800.
     rme(9,is) = 1700.
     rme(10,is) = 1680.
     print(sname(is)+" "+rme(:,is))
  end if

  end do 


  print("--------- valid ann sites -----------") 
  do is = 0,ns-1 
    ;print("-") 
    ;print("obs: "+sname(is)+"  "+rme(:,is)+"   "+rma(:,is)+"   "+rmd(:,is)) 
    ;print("-") 
    if(.not.ismissing(sme(is)) .and. .not. ismissing(rme(0,is)) ) then 
       if(abs(sme(is)-rme(0,is)).le.1.) then 
          rme(:,is) = -999 
          rma(:,is) = -999 
          rmb(:,is) = -999 
          rmc(:,is) = -999 
          rmd(:,is) = -999 
          knot(is)  = 0
          if(kbad(is).eq.1) then  
             print(strnice(sname(is),20)+"") 
          end if 
          kbad(is) = 0 
       end if 
    end if 
  end do 
  print("") 
  print("") 

  print("--------- valid sites -----------")
  do is = 0,ns-1 
    if(knot(is).ne.0) then 
       print(strnice(sname(is),20)+"       "+sprintf("%9.3f",sme(is))+"    "+sprintf("%9.3f",smd(is)) ) 
    else
       kbad(is) = 0 
    end if 
  end do 
  print("") 
  print("") 
  print("number of valid sites: "+sprinti("%6.0i",sum(kbad)) ) 
  print("") 
  print("") 
  print("------------------------------------------------------------------") 

;;----------------------------------------------------------------------
;; recalculate annual mean 
;;----------------------------------------------------------------------
  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  
 
;;----------------------------------------------------------------------
;; mask for DJF, JJA values 
;;----------------------------------------------------------------------
  txa = rma 
  tya = rma
  tza = rma
  toa = rma ;;var_onlySON(rma,0) 

  do is = 0,ns-1 

  if(lat(is).gt.0.) then ;; NH 
     txa(:,is) = var_onlyDJF(rma(:,is),0)
     tya(:,is) = var_onlyJJA(rma(:,is),0) 
     tza(:,is) = var_onlyMAM(rma(:,is),0) 
     toa(:,is) = var_onlySON(rma(:,is),0) 
  else
     txa(:,is) = var_onlyJJA(rma(:,is),0)
     tya(:,is) = var_onlyDJF(rma(:,is),0) 
     tza(:,is) = var_onlySON(rma(:,is),0) 
     toa(:,is) = var_onlyMAM(rma(:,is),0) 
  end if 

  end do

  txb = rmc 
  tyb = rmc
  tzb = rmc
  tob = rmc ;;var_onlySON(rmc,0)

  do is = 0,ns-1

  if(lat(is).gt.0.) then ;; NH 
     txb(:,is) = var_onlyDJF(rmc(:,is),0)
     tyb(:,is) = var_onlyJJA(rmc(:,is),0)
     tzb(:,is) = var_onlyMAM(rmc(:,is),0)
     tob(:,is) = var_onlySON(rmc(:,is),0)
  else
     txb(:,is) = var_onlyJJA(rmc(:,is),0)
     tyb(:,is) = var_onlyDJF(rmc(:,is),0)
     tzb(:,is) = var_onlySON(rmc(:,is),0)
     tob(:,is) = var_onlyMAM(rmc(:,is),0)
  end if

  end do
 
  txc = rmd
  tyc = rmd
  tzc = rmd
  toc = rmd ;;var_onlySON(rmd,0)

  do is = 0,ns-1

  if(lat(is).gt.0.) then ;; NH 
     txc(:,is) = var_onlyDJF(rmd(:,is),0)
     tyc(:,is) = var_onlyJJA(rmd(:,is),0)
     tzc(:,is) = var_onlyMAM(rmd(:,is),0)
     toc(:,is) = var_onlySON(rmd(:,is),0)
  else
     txc(:,is) = var_onlyJJA(rmd(:,is),0)
     tyc(:,is) = var_onlyDJF(rmd(:,is),0)
     tzc(:,is) = var_onlySON(rmd(:,is),0)
     toc(:,is) = var_onlyMAM(rmd(:,is),0)
  end if

  end do
 
  txe = rme 
  tye = rme
  tze = rme
  toe = rme ;;var_onlySON(rme,0) 

  do is = 0,ns-1

  if(lat(is).gt.0.) then ;; NH 
     txe(:,is) = var_onlyDJF(rme(:,is),0)
     tye(:,is) = var_onlyJJA(rme(:,is),0)
     tze(:,is) = var_onlyMAM(rme(:,is),0)
     toe(:,is) = var_onlySON(rme(:,is),0)
  else
     txe(:,is) = var_onlyJJA(rme(:,is),0)
     tye(:,is) = var_onlyDJF(rme(:,is),0)
     tze(:,is) = var_onlySON(rme(:,is),0)
     toe(:,is) = var_onlyMAM(rme(:,is),0)
  end if

  end do

;;----------------------------------------------------------------------
;; arrage variables  
;;----------------------------------------------------------------------

  rmea = rme 
  rmeb = rme 
  rmec = rme 
  rmed = rme 

  rmaa = rma 
  rmab = rma 
  rmac = rma 
  rmad = rma 

  rmca = rmc 
  rmcb = rmc 
  rmcc = rmc 
  rmcd = rmc 

  rmda = rmd 
  rmdb = rmd
  rmdc = rmd 
  rmdd = rmd 

  do is = 0,ns-1 
  if(ieur(is).ne.0) then 
     rmea(:,is) = -999 
     rmaa(:,is) = -999 
     rmca(:,is) = -999 
     rmda(:,is) = -999 
  end if  
  print(sname(is)+"   "+max(rme(:,is))) 
  end do 

  do is = 0,ns-1
  if(ichi(is).ne.0) then
     rmeb(:,is) = -999
     rmab(:,is) = -999
     rmcb(:,is) = -999
     rmdb(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  do is = 0,ns-1
  if(iusa(is).ne.0) then
     rmec(:,is) = -999
     rmac(:,is) = -999
     rmcc(:,is) = -999
     rmdc(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  do is = 0,ns-1
  if(isam(is).ne.0) then
     rmed(:,is) = -999
     rmad(:,is) = -999
     rmcd(:,is) = -999
     rmdd(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  x1  = ndtooned(rme(:,:)) 
  y1  = ndtooned(rma(:,:)) 
  y2  = ndtooned(rmc(:,:)) 
  y3  = ndtooned(rmd(:,:)) ;; a shift !!!  

  x1a  = ndtooned(rmea(:,:)) 
  y1a  = ndtooned(rmaa(:,:)) 
  y2a  = ndtooned(rmca(:,:)) 
  y3a  = ndtooned(rmda(:,:)) ;; a shift !!!  

  x1b  = ndtooned(rmeb(:,:)) 
  y1b  = ndtooned(rmab(:,:)) 
  y2b  = ndtooned(rmcb(:,:)) 
  y3b  = ndtooned(rmdb(:,:)) ;; a shift !!!  

  x1c  = ndtooned(rmec(:,:)) 
  y1c  = ndtooned(rmac(:,:)) 
  y2c  = ndtooned(rmcc(:,:)) 
  y3c  = ndtooned(rmdc(:,:)) ;; a shift !!!  

  x1d  = ndtooned(rmed(:,:)) 
  y1d  = ndtooned(rmad(:,:)) 
  y2d  = ndtooned(rmcd(:,:)) 
  y3d  = ndtooned(rmdd(:,:)) ;; a shift !!!  


  ;; seasons 

  xa1 = ndtooned(txe(:,:)) ;; OBS ;; DJF 
  ya1 = ndtooned(txa(:,:)) ;; WCRP1995
  ya2 = ndtooned(txb(:,:)) ;; SH1998 
  ya3 = ndtooned(txc(:,:)) ;; MERGE

  xb1 = ndtooned(tye(:,:)) ;; JJA
  yb1 = ndtooned(tya(:,:))
  yb2 = ndtooned(tyb(:,:))
  yb3 = ndtooned(tyc(:,:))

  xc1 = ndtooned(tze(:,:)) ;; MAM
  yc1 = ndtooned(tza(:,:))
  yc2 = ndtooned(tzb(:,:))
  yc3 = ndtooned(tzc(:,:))

  xd1 = ndtooned(toe(:,:)) ;; SON
  yd1 = ndtooned(toa(:,:))
  yd2 = ndtooned(tob(:,:))
  yd3 = ndtooned(toc(:,:))

;;----------------------------------------------------------------------
;; mask negative values again 
;;----------------------------------------------------------------------
  x1  = var_maskneg(x1,0) 
  y1  = var_maskneg(y1,0) 
  y2  = var_maskneg(y2,0) 
  y3  = var_maskneg(y3,0) 

  x1a  = var_maskneg(x1a,0) 
  y1a  = var_maskneg(y1a,0) 
  y2a  = var_maskneg(y2a,0) 
  y3a  = var_maskneg(y3a,0) 

  x1b  = var_maskneg(x1b,0) 
  y1b  = var_maskneg(y1b,0) 
  y2b  = var_maskneg(y2b,0) 
  y3b  = var_maskneg(y3b,0) 

  x1c  = var_maskneg(x1c,0) 
  y1c  = var_maskneg(y1c,0) 
  y2c  = var_maskneg(y2c,0) 
  y3c  = var_maskneg(y3c,0) 

  x1d  = var_maskneg(x1d,0) 
  y1d  = var_maskneg(y1d,0) 
  y2d  = var_maskneg(y2d,0) 
  y3d  = var_maskneg(y3d,0) 

  xa1 = var_maskneg(xa1,0) 
  ya1 = var_maskneg(ya1,0) 
  ya2 = var_maskneg(ya2,0) 
  ya3 = var_maskneg(ya3,0) 

  xb1 = var_maskneg(xb1,0) 
  yb1 = var_maskneg(yb1,0) 
  yb2 = var_maskneg(yb2,0) 
  yb3 = var_maskneg(yb3,0) 

  xc1 = var_maskneg(xc1,0) 
  yc1 = var_maskneg(yc1,0) 
  yc2 = var_maskneg(yc2,0) 
  yc3 = var_maskneg(yc3,0) 

  xd1 = var_maskneg(xd1,0) 
  yd1 = var_maskneg(yd1,0) 
  yd2 = var_maskneg(yd2,0) 
  yd3 = var_maskneg(yd3,0) 
  ;print("x1,y1,y2: "+x1+"   "+y1+"   "+y2) 

  fy1 = var_p2(x1,y1,0) 
  fy2 = var_p2(x1,y2,0) 
  fy3 = var_p2(x1,y3,0) 
  
  fy1a = var_p2(x1a,y1a,0) 
  fy2a = var_p2(x1a,y2a,0) 
  fy3a = var_p2(x1a,y3a,0) 
  
  fy1b = var_p2(x1b,y1b,0) 
  fy2b = var_p2(x1b,y2b,0) 
  fy3b = var_p2(x1b,y3b,0) 
  
  fy1c = var_p2(x1c,y1c,0) 
  fy2c = var_p2(x1c,y2c,0) 
  fy3c = var_p2(x1c,y3c,0) 
  
  fy1d = var_p2(x1d,y1d,0) 
  fy2d = var_p2(x1d,y2d,0) 
  fy3d = var_p2(x1d,y3d,0) 
  
  fya1 = var_p2(xa1,ya1,0) 
  fya2 = var_p2(xa1,ya2,0) 
  fya3 = var_p2(xa1,ya3,0) 

  fyb1 = var_p2(xb1,yb1,0) 
  fyb2 = var_p2(xb1,yb2,0) 
  fyb3 = var_p2(xb1,yb3,0) 

  fyc1 = var_p2(xc1,yc1,0) 
  fyc2 = var_p2(xc1,yc2,0) 
  fyc3 = var_p2(xc1,yc3,0) 

  fyd1 = var_p2(xd1,yd1,0) 
  fyd2 = var_p2(xd1,yd2,0) 
  fyd3 = var_p2(xd1,yd3,0) 

  print(""+fy1+" "+fy2+" "+fy3)
  print(""+fya1+" "+fya2+" "+fya3)
  print(""+fyb1+" "+fyb2+" "+fyb3)
  print(""+fyc1+" "+fyc2+" "+fyc3)
  print(""+fyd1+" "+fyd2+" "+fyd3)
 
;;----------------------------------------------------------------------
;; plotvar  
;;----------------------------------------------------------------------
  yy = new((/3,dimsizes(y1)/),"float") 

  yy(0,:) = y1(:) 
  yy(1,:) = y2(:) 
  yy(2,:) = y3(:) 

  ya = new((/3,dimsizes(y1)/),"float") 
  ya(0,:) = ya1(:) 
  ya(1,:) = ya2(:) 
  ya(2,:) = ya3(:) 

  yb = new((/3,dimsizes(y1)/),"float") 
  yb(0,:) = yb1(:) 
  yb(1,:) = yb2(:) 
  yb(2,:) = yb3(:) 

  yc = new((/3,dimsizes(y1)/),"float") 
  yc(0,:) = yc1(:) 
  yc(1,:) = yc2(:) 
  yc(2,:) = yc3(:) 

  yd = new((/3,dimsizes(y1)/),"float") 
  yd(0,:) = yd1(:) 
  yd(1,:) = yd2(:) 
  yd(2,:) = yd3(:) 


;;----------------------------------------------------------------------
;; correlation 
;;----------------------------------------------------------------------
  fr1 = esccr(x1,y1,0)
  fr2 = esccr(x1,y2,0)
  fr3 = esccr(x1,y3,0)
 
  fra1 = esccr(xa1,ya1,0)
  fra2 = esccr(xa1,ya2,0)
  fra3 = esccr(xa1,ya3,0)
 
  frb1 = esccr(xd1,yd1,0)
  frb2 = esccr(xd1,yd2,0)
  frb3 = esccr(xd1,yd3,0)
 
  frc1 = esccr(xc1,yc1,0)
  frc2 = esccr(xc1,yc2,0)
  frc3 = esccr(xc1,yc3,0)
 
  frd1 = esccr(xb1,yb1,0)
  frd2 = esccr(xb1,yb2,0)
  frd3 = esccr(xb1,yb3,0)
 
  print("") 
  print(""+fr1+" "+fr2+" "+fr3)
  print(""+fra1+" "+fra2+" "+fra3)
  print(""+frb1+" "+frb2+" "+frb3)
  print(""+frc1+" "+frc2+" "+frc3)
  print(""+frd1+" "+frd2+" "+frd3)

  cr = new(np,"float") 
  p2 = new(np,"float") 
  cr = 0. 
  p2 = 0. 

  cr(0) = fr1 
  cr(1) = fr2 
  cr(2) = fr3 
  cr(3) = fra1 
  cr(4) = fra2 
  cr(5) = fra3 
  cr(6) = frb1 
  cr(7) = frb2 
  cr(8) = frb3 

  p2(0) = fy1
  p2(1) = fy2
  p2(2) = fy3
  p2(3) = fya1
  p2(4) = fya2
  p2(5) = fya3
  p2(6) = fyb1
  p2(7) = fyb2
  p2(8) = fyb3

;;----------------------------------------------------------------------
;; open wks and plot 
;;----------------------------------------------------------------------
 

  wks = gsn_open_wks(plotform,plotname)
  gsn_define_colormap(wks,colormap)


;;----------------------------------------------------------------------
;; plot setting 
;;----------------------------------------------------------------------
  res          = True
  res@gsnDraw  = False
  res@gsnFrame = False

  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@tiXAxisFontHeightF     = axisfont 
  res@tiYAxisFontHeightF     = axisfont 
  res@tmYLLabelFontHeightF   = axisfont
  res@tmXBLabelFontHeightF   = axisfont

  res@gsnLeftString              = LeftString
  res@gsnRightString             = RightString 
  res@gsnLeftStringFontHeightF   = axisfont 
  res@gsnRightStringFontHeightF  = axisfont 
 
 
  ;;res@mpShapeMode  = "FreeAspect"
  ;;res@vpWidthF      = 0.6
  ;;res@vpHeightF     = 0.6

  
  if(L_log.eq.0) 
    res@trXLog     = True
    res@trYLog     = True
    res@trXMinF    = 1
    res@trXMaxF    = 50000
    res@trYMinF    = 1
    res@trYMaxF    = 50000
  else
    res@trXLog     = False
    res@trYLog     = False
    res@trXMinF    = 0.
    res@trXMaxF    = 5000
    res@trYMinF    = 0.
    res@trYMaxF    = 5000
  end if
  
  res@trXReverse     = False

;;----------------------------------------------------------------------
;; marker setting 
;;----------------------------------------------------------------------
  res@xyMarkLineModes   = "Markers"
  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize
  res@xyMarkerColor     = icla ;;"blue"

;;----------------------------------------------------------------------
;; raw plot 
;;----------------------------------------------------------------------

  ip = 0 

  res@tiXAxisString       = "" 
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "ALL" 
  plot(ip) = gsn_csm_xy(wks,x1,yy(0,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiYAxisString       = "" 
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "ALL" 
  plot(ip) = gsn_csm_xy(wks,x1,yy(1,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiYAxisString       = "" 
  res@gsnLeftString       = abcd(ip) + ") MERGED" 
  res@gsnRightString      = "ALL" 
  plot(ip) = gsn_csm_xy(wks,x1,yy(2,:),res)

  res@xyMarkerColor     = iclb ;;"green"
;;; 

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "DJF" 
  plot(ip) = gsn_csm_xy(wks,xa1,ya(0,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "DJF" 
  plot(ip) = gsn_csm_xy(wks,xa1,ya(1,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") MERGED" 
  res@gsnRightString      = "DJF" 
  plot(ip) = gsn_csm_xy(wks,xa1,ya(2,:),res)

  res@xyMarkerColor     = iclc ;;"brown"
;;; 

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "JJA" 
  plot(ip) = gsn_csm_xy(wks,xb1,yb(0,:),res)

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "JJA" 
  plot(ip) = gsn_csm_xy(wks,xb1,yb(1,:),res)

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") MERGED" 
  res@gsnRightString      = "JJA" 
  plot(ip) = gsn_csm_xy(wks,xb1,yb(2,:),res)


;;----------------------------------------------------------------------
;; print P2 and CR values on plot 
;;----------------------------------------------------------------------
  txres               = True
  txres@txFontHeightF = textfont 
  txres@txJust        = "CenterLeft" 

  texta = new(np,graphic)  
  textb = new(np,graphic)  

  do ip = 0,np-1
     texta(ip) = gsn_add_text(wks,plot(ip),"P2: "+sprintf("%3.2f",p2(ip))+"%",1000.,18.,txres)
     textb(ip) = gsn_add_text(wks,plot(ip),"CR: "+sprintf("%3.2f",cr(ip)),1000.,10.,txres)
  end do 

;;----------------------------------------------------------------------
;; polygon 
;;----------------------------------------------------------------------
  xp    = new(6,"float")
  yp    = new(6,"float")

  xp(0) = 1.0 * res@trXMinF  
  yp(0) = 1.0 * res@trYMinF  
  xp(1) = 2.0 * res@trXMinF 
  yp(1) = 1.0 * res@trYMinF  
  xp(2) = 1.0 * res@trXMaxF 
  yp(2) = 0.5 * res@trYMaxF 
  xp(3) = 1.0 * res@trXMaxF 
  yp(3) = 1.0 * res@trYMaxF 
  xp(4) = 0.5 * res@trXMaxF  
  yp(4) = 1.0 * res@trYMaxF  
  xp(5) = 1.0 * res@trXMinF  
  yp(5) = 2.0 * res@trYMinF  

 ;gsres                   = True                        ; poly res
 ;gsres@tfPolyDrawOrder   = "Predraw"                   ; draw this first
 ;gsres@gsFillColor       = "grey"                 ; color chosen
 ;duma = gsn_add_polygon(wks,plot0,xp,yp,gsres)

;;----------------------------------------------------------------------
;; polyline 
;;----------------------------------------------------------------------
  resl = True
  resl@gsLineColor = "black"
  resl@gsLineLabelString= ""

  xpts = (/1.,1./)
  ypts = (/1.,1./)

  ypts(0) = res@trYMinF
  ypts(1) = res@trYMaxF
  xpts(0) = res@trXMinF
  xpts(1) = res@trXMaxF

  resl@gsLineThicknessF = 0.1
  resl@gsLineDashPattern = 0

  dunaa=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  dunab=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  dunac=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  dunad=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunae=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  dunaf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  dunag=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunah=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunai=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)

  resl@gsLineThicknessF = 1.0
  resl@gsLineDashPattern = 11

  ypts(0) = res@trYMinF
  ypts(1) = 0.5 * res@trYMaxF
  xpts(0) = 2. * res@trXMinF
  xpts(1) = res@trXMaxF

  dunba=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  dunbb=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  dunbc=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  dunbd=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunbe=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  dunbf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  dunbg=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunbh=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunbi=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)

  ypts(0) = 2. * res@trYMinF
  ypts(1) = res@trYMaxF
  xpts(0) = res@trXMinF
  xpts(1) = 0.5 * res@trXMaxF

  dunca=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  duncb=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  duncc=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  duncd=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunce=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  duncf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  duncg=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunch=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunci=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)
  ;duncj=gsn_add_polyline(wks,plot(9),xpts(0:1),ypts(0:1),resl)
  ;dunck=gsn_add_polyline(wks,plot(10),xpts(0:1),ypts(0:1),resl)
  ;duncl=gsn_add_polyline(wks,plot(11),xpts(0:1),ypts(0:1),resl)
  ;duncm=gsn_add_polyline(wks,plot(12),xpts(0:1),ypts(0:1),resl)
  ;duncn=gsn_add_polyline(wks,plot(13),xpts(0:1),ypts(0:1),resl)
  ;dunco=gsn_add_polyline(wks,plot(14),xpts(0:1),ypts(0:1),resl)

;;----------------------------------------------------------------------
;; pannel 
;;----------------------------------------------------------------------
  ResP                            = True
;;ResP@gsnMaximize                = True
;;ResP@gsnPaperMargin             = 3 
;;ResP@gsnPanelXWhiteSpacePercent = 1
;;ResP@gsnPanelYWhiteSpacePercent = 1 

  gsn_panel (wks,plot, (/3,3/),ResP)


;;----------------------------------------------------------------------
;; output bad site index (0) 
;;----------------------------------------------------------------------
  asciiwrite("kbad.dat",kbad) 

  system("ps2pdf fig_scatter_color.eps")  
  ;draw(plot)     
  ;frame(wks)


